## This file creates all figures in the main body and supplementary materials.
setwd("/Users/agunderson/Dropbox/Book Project")

#clear working environment
rm(list=ls())

library(ggplot2)
library(readxl)
library(tidyverse)
library(geofacet)
library(reshape2)

##################################################################
# Figure 1 in the paper --------------
##################################################################  
lawsuits <- read_csv("./Data/fjcIVtermandfile.csv") %>% 
  filter(Year>1984)%>%
  group_by(Year,Category) %>% 
  dplyr::summarise(sumlawsuitsall = sum(sumlawsuitsall))%>%
  as.data.frame() 

ggplot(lawsuits,aes(y=sumlawsuitsall,x=Year,colour=Category)) + geom_point()+
  xlab("Year")+ylab("Number of Lawsuits")+
  ggtitle("Prisoner Lawsuits All States, 1986 to 2016")+
  theme(plot.title = element_text(hjust = 0.5)) + theme_bw()+scale_colour_grey()
ggsave("./Figures/Lawsuitsovertime.pdf")

##################################################################
# Figure 2 in the paper --------------
##################################################################  
incarrate <- read_csv("./Data/nationalincarrate.csv") 
incarrate.melt <- gather(incarrate, Jurisdiction,incarrate, Federal:Private, factor_key=TRUE)
ggplot(incarrate.melt,aes(x=year,y=incarrate,colour=Jurisdiction,shape=Jurisdiction))+
  geom_point()+theme_bw()+xlab("Year")+
  ylab("Incarceration Rate \n (Per 100,000 Resident Population)")+xlim(1984,2012) +
  scale_colour_grey(start = 0, end = .9)
ggsave("./Figures/incarrate.pdf")

##################################################################
# Figure 3 in the paper --------------
##################################################################  
fjcandprivatization <- read_csv("./Data/POPdata.csv")

ggplot(fjcandprivatization, aes(year, logDC)) +
  geom_line() +
  facet_geo(~ jurisdiction, grid = "us_state_grid3", label = "code") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  labs(title = "Logged Number of Private Prison Inmates, 1986 - 2016",
       x = "Year",
       y = "Logged Number of Private Inmates") +
  theme(strip.text.x = element_text(size = 6,face = "bold"))+theme_bw()
ggsave("./Figures/maplogppinmates.pdf",height=8,width = 12)

##################################################################
# Supplementary Figure 1 --------------
##################################################################  
fjcIV <- read_csv("./Data/fjclength.csv")

ggplot(fjcIV, aes(x=filetotermdiff.year)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  geom_density(alpha=.2, fill="darkgray") +theme_bw() +xlab("Year to Termination")+ylab("Density")
ggsave("./Figures/lengthtotermination.pdf")

##################################################################
# Supplementary Figures 2-4 --------------
##################################################################  
ggplot(fjcandprivatization, aes(year, propDC)) +
  geom_line() +
  facet_geo(~ jurisdiction, grid = "us_state_grid3", label = "code") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  labs(title = "Proportion of Private Prison Inmates, 1986 - 2016",
       x = "Year",
       y = "Prop. of Private Inmates") +
  theme(strip.text.x = element_text(size = 6,face = "bold"))+theme_bw()
ggsave("./Figures/mappropppinmates.pdf",height=8,width = 12)

ggplot(fjcandprivatization, aes(year, sumfacilities)) +
  geom_line() +
  facet_geo(~ jurisdiction, grid = "us_state_grid3", label = "code") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  labs(title = "Number of Private Facilities, 1986 - 2016",
       x = "Year",
       y = "Number of Private Facilities") +
  theme(strip.text.x = element_text(size = 6,face = "bold"))+theme_bw()
ggsave("./Figures/mapppfacilities.pdf",height=8,width = 12)

fjcandprivatization <- fjcandprivatization %>%
  mutate(total.prison.pop = (incar.rate.y*total.pop)/100000,
    lawsuitsperprisoner = (sumlawsuitsall / total.prison.pop))
ggplot(fjcandprivatization, aes(year, lawsuitsperprisoner)) +
  geom_line() +
  facet_geo(~ jurisdiction, grid = "us_state_grid3", label = "code") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  labs(title = "Lawsuits Per Prisoner, 1986 - 2016",
       x = "Year",
       y = "Lawsuits per Prisoner") +
  theme(strip.text.x = element_text(size = 6,face = "bold"))+theme_bw()
ggsave("./Figures/maplawsuitsperprisoner.pdf",height=8,width = 12)

##################################################################
# Supplementary Figures 5-6 --------------
################################################################## 
latlongfacilities <- read_csv("./raw-data/latlongfacilities.csv") %>% 
  nest(MIN_year, MAX_year) %>%
  mutate(data = map(data, ~seq(unique(.x$MIN_year), unique(.x$MAX_year), 1))) %>%
  unnest(data) %>% rename(year=data)%>%
  dplyr::filter(., !grepl("Puerto Rico",DisCode))

juris.map <- latlongfacilities %>%
  group_by(FIRST_Location,POINT_X,POINT_Y,year)%>%
  summarise(sumfacilities=n())  

map.function <- function(datayear){
  usa <- map_data("state")
  ditch_the_axes <- theme(axis.text = element_blank(),axis.line = element_blank(),axis.ticks = element_blank(),
                          panel.border = element_blank(),panel.grid = element_blank(),axis.title = element_blank())
  citiesmap <- ggplot() + geom_polygon(data = usa, aes(x=long, y = lat,group=group),color="darkgray",fill="NA")
  citiesmap+geom_point(data = juris.map[which(juris.map$year==datayear),], aes(x = POINT_X, y = POINT_Y,
                                                                               size=sumfacilities),alpha=0.7,shape = 21, colour = "black", fill = "gray70", 
                       stroke = 0.4)+
    theme_bw()+ditch_the_axes+ theme(legend.position="none")  
  ggsave(paste0("./Figures/",datayear,".pdf"),width=11,height=7)
}

map.function("1986")
map.function("1996")
map.function("2006")
map.function("2016") 

##################################################################
# Supplementary Figure 7 --------------
################################################################## 
customers <- fjcandprivatization %>%
  select(year,sumfacilitiesonecustomer,sumfacilitiesmorethanonecustomer) %>%
  replace(is.na(.), 0) %>%
  group_by(year) %>%
  summarise(sumfacilitiesonecustomerall = sum(sumfacilitiesonecustomer),
            sumfacilitiesmorethanonecustomerall = sum(sumfacilitiesmorethanonecustomer)) %>% 
  rename(Year=year,`One Customer`= sumfacilitiesonecustomerall, 
         `More Than One Customer`= sumfacilitiesmorethanonecustomerall) %>%
  melt(.,id="Year") %>%
  rename(Customers=variable) 

ggplot(customers,aes(x=Year,y=value,color=Customers)) + geom_point()+ 
  stat_smooth(aes(group=Customers)) + ylab("Number of Facilities") + 
  ggtitle("Distribution of Number of Customers \n by Year") + theme_bw() + scale_color_grey(start=.4,end=.7)
ggsave("./Figures/numcustomers.pdf")

##################################################################
# Supplementary Figure 8 --------------
##################################################################  
bjsandSEC <- fjcandprivatization %>%
  select(year,private.pop,designcapacity2)%>%
  replace(is.na(.), 0) %>%
  group_by(year) %>%
  summarise(BJS = sum(private.pop),
            "Original Data" = sum(designcapacity2)) 

bjsandSEClong <- melt(bjsandSEC, id="year") %>% 
  rename(Source=variable)

cor(bjsandSEC$BJS,bjsandSEC$`Original Data`) #0.8436206

ggplot(bjsandSEClong,
       aes(x=year, y=value, colour=Source)) + geom_point()+ 
  stat_smooth(aes(group=Source),method="lm") + xlim(1999,2015) + ylim(45000,105000)+
  ylab("Sum of Private Prisoners") + xlab("Year") + 
  ggtitle("Correlation Between BJS Data \n and Original Dataset") + theme_bw()+ scale_color_grey(start=.4,end=.7)
ggsave("./Figures/bjsdatacorrelation.pdf")


